# Replication for "Terrorism and the Rise of Right-Wing Content in Israeli Books"
# Tamar Mitts (2017)

#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#

# This code replicates the figures in Mitts (2017). 
# The main analysis is done in Stata. See relevant Stata code in the Replication folder.

#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#


rm(list=ls())

library(foreign)
library(wordcloud)
library(gplots)

setwd("~/Downloads/Replication/")

# Figure 1: Rising support for right-wing parties over time

ines_data = read.csv("Data/bloc_vote_for.csv", as.is=T)
ines_data_t = t(ines_data[,2:4])

pdf("Figures/ines_over_time_barplot.pdf", width=8, height=5)
par(mar=c(3,3,2,2))
barplot(height= as.matrix(ines_data_t), beside=T, border="black", col="black", names.arg= as.character(ines_data$year), ylim=c(0,1), density=c(20,35,100), angle=c(135, 45, 45))
X = c(9,20)
Y = c(0,1)
rect(X[1], Y[1], X[2], Y[2], border = "gray80", col = "gray80")
barplot(height= as.matrix(ines_data_t), beside=T, axes=F, xlab="", ylab="", border="black", col="black", names.arg=as.character(ines_data$year), add=T, ylim=c(0,1), xaxt="n", density=c(20, 35, 100), angle=c(135, 45, 45))
text(14.5, 0.9, labels="2nd \nIntifada", cex=1)
legend(28.5, 1, legend=c("Right", "Center", "Left"), cex=0.8, bg="white", box.col="white", density=c(100, 35, 20), angle=c(45, 45, 135))
dev.off()

# Figure 2: Decrease in public support for negotiations over time

peace_index = read.csv("Data/Negotiation_index.csv", as.is=T)

neg_index_by_year <- aggregate(x= peace_index$Negotiations.Index, by= list(peace_index$Year), FUN=mean, na.rm=T)

pdf(file="Figures/negotiation_index.pdf", width=8, height=5)
par(mar=c(4,4,2,2))
plot(neg_index_by_year$Group.1, neg_index_by_year$x, type="l", ylab="Negotitations Index", xlab="Year", col="white")
rect(xleft=2000, ybottom=44.5, xright=2005, ytop=64.5, col="lightgray", border=NA)
lines(neg_index_by_year$Group.1, neg_index_by_year$x, col="red", lwd=2)
text(2002.5, 62, labels="Second Intifada")
abline(h=50)
dev.off()

# Figure 3: Phrases most strongly linked the ideology of political blocs

data.right <- read.csv("Data/top100-right-phrase.csv", header=TRUE, as.is=T)
right.heb <- readLines("Data/right-phrases_heb.txt")
right.eng <- readLines("Data/right-phrases_eng.txt")
data.right <- data.right[data.right$top.right.phrases %in% right.heb,]
data.right$right.eng <- right.eng
data.right <- data.right[order(data.right$freq, decreasing=T),]
data.right <- data.right[1:50,]

data.center <- read.csv("Data/top100-center-phrase.csv", header=TRUE, as.is=T)
center.heb <- readLines("Data/center-phrases_heb.txt")
center.eng <- readLines("Data/center-phrases_eng.txt")
data.center <- data.center[data.center$top.center.phrases %in% center.heb,]
data.center$center.eng <- center.eng
data.center <- data.center[order(data.center$freq, decreasing=T),]
data.center <- data.center[1:50,]

data.left <- read.csv("Data/top100-left-phrase.csv", header=TRUE, as.is=T)
left.heb <- readLines("Data/left-phrases_heb.txt")
left.eng <- readLines("Data/left-phrases_eng.txt")
data.left <- data.left[data.left$top.left.phrases %in% left.heb,]
data.left$left.eng <- left.eng
data.left <- data.left[order(data.left$freq, decreasing=T),]
data.left <- data.left[1:50,]

pdf("Figures/word_cloud_left.pdf", width=5, height=5)
size.phrases <- 1.7
par(mar=c(1,1,1,1))
wordcloud(words= data.left$left.eng, freq= data.left$freq, scale=c(1.8,1), rot.per=0, random.order=F, colors=colorRampPalette(c('blue','black'))(length(data.left$top.left.phrases)), min.freq=0)
text(0.5,1, labels="Left-wing", font=2, cex= size.phrases)
dev.off()

pdf("Figures/word_cloud_center.pdf", width=5, height=5)
size.phrases <- 1.7
par(mar=c(1,1,1,1))
wordcloud(words= data.center$center.eng, freq= data.center$freq, scale=c(1.8,1), rot.per=0, random.order=F, colors=colorRampPalette(c('purple','black'))(length(data.center$top.center.phrases)), min.freq=0)
text(0.5,1, labels="Center", font=2, cex= size.phrases)
dev.off()

pdf("Figures/word_cloud_right.pdf", width=5, height=5)
size.phrases <- 1.7
par(mar=c(1,1,1,1))
wordcloud(words= data.right$right.eng , freq= data.right$freq, scale=c(1.8,1), rot.per=0, random.order=F, colors=colorRampPalette(c('brown','black'))(length(data.right$top.right.phrases)), min.freq=0)
text(0.5,1, labels="Right-wing", font=2, cex= size.phrases)
dev.off()

# Figure 5: Summary statistics: Phrase frequency over time

sumstats_table = read.csv("Data/summary_statistics.csv", as.is=T)
sumstats_table$freqX10000 = sumstats_table$Frequency_mean*10000

pdf(file="Figures/sumstats_figure.pdf", width=9, height=4)
par(mar=c(3,5,2,5))
with(sumstats_table, barplot(Phrases, col="gray90",names.arg= sumstats_table$Year, cex.names =0.8, border=NA, cex.axis=0.8, ylab="Number of phrases"))
par(new = T)
with(sumstats_table, scatter.smooth(Year, freqX10000, col=NA, bty="n", ylab=NA, xlab=NA, xaxt="n", yaxt="n"))
lines(loess.smooth(sumstats_table$Year, sumstats_table$freqX10000), col="gray20", lwd=3, cex.axis=0.8)
axis(side = 4, cex.axis=0.8)
mtext(side = 4, line = 3, "Average frequency (x 10,000)")
legend("topleft", legend=c("Number of phrases", "Average frequency"), pch = c(15, NA) ,pt.cex=c(2,NA),  lty = c(NA, 1), lwd = c(NA, 3), col = c("gray90", "black"), bg="white", box.col="white")
dev.off()

# Figure 6: Israeli casualties in the Intifadas

israelis.killed = read.csv("Data/Israeli_casualties.csv", as.is=T)

pdf(file="Figures/israeli_casualties_figure.pdf", width=8, height=4)
par(mar=c(5,5,1,1))
plot(israelis.killed$Year, israelis.killed$Total.killed, type="l", main="", xlab="Year", ylab="Number of casualties", col=NA, frame.plot=F, xlim=c(1985,2010), cex.axis=0.8)
rect(xleft=1987, ybottom=0, xright=1991, ytop=max(israelis.killed$Total.killed)+10, col="gray80", border=NA)
rect(xleft=2000, ybottom=0, xright=2005, ytop=max(israelis.killed$Total.killed)+10, col="gray80", border=NA)
lines(israelis.killed$Year, israelis.killed$Total.killed)
text(1989, 50, labels="First Intifada", cex=0.8)
text(2002.5, 30, labels="Second Intifada", cex=0.8)
dev.off()

# Figure 7: Terrorism and right-wing content in books

coeffs.results <- read.csv("Data/main_results_coefs_for_graph.csv")
pct.change <- read.csv("Data/predicted_values_results.csv")
coeffs.results_right <- coeffs.results[coeffs.results$bloc=="right",]
results_right_second_intif_peak <- coeffs.results_right[1:3,]
results_right_second_intif_peak$pct_change <- pct.change$pct_change_right[1:3]
results_right_second_intif_casualties <- coeffs.results_right[4:6,]
results_right_second_intif_casualties$pct_change <- pct.change$pct_change_right[4:6]

pdf(file="Figures/main_results_peak_secondint_right.pdf", width=9, height=5)
par(mar=c(5,4,3,4))
labels.peak <- c("[-6,-1], \n[1,6]", "[-6,-4], \n[4,6]", "[-3,-1], \n[1,3]")
labels.casualties <- c("Cumulative", "MA5", "L5")
par(mfrow=c(1,2))
plotCI(x= results_right_second_intif_peak$estimate, ui= results_right_second_intif_peak$max95, li= results_right_second_intif_peak$min95, type="p", pt.bg= "red", col="red", pch=21, gap=0, xaxt="n", ylab="Coefficient", xlab="", cex.main=1, cex.axis=0.9, cex=1, sfrac=0, ylim=c(-1,1), xlim=c(0.5,3.5), main="Peak of Terrorism")
abline(h=0, col="black")
axis(side=1, at=1:3, labels= labels.peak, cex.axis=0.8, padj=1)
text(1.3,results_right_second_intif_peak$estimate[1], labels=paste0(round(results_right_second_intif_peak$pct_change[1], 2),"%"), col="red", cex=0.8)
text(2.3,results_right_second_intif_peak$estimate[2], labels=paste0(round(results_right_second_intif_peak$pct_change[2], 2),"%"), col="red", cex=0.8)
text(3.3,results_right_second_intif_peak$estimate[3], labels=paste0(round(results_right_second_intif_peak$pct_change[3], 2),"%"), col="gray60", cex=0.8)

plotCI(x= results_right_second_intif_casualties$estimate, ui= results_right_second_intif_casualties$max95, li= results_right_second_intif_casualties$min95, type="p", pt.bg= "red", col="red", pch=21, gap=0, ylab="Coefficient", xaxt="n", xlab="", cex.main=1, cex.axis=0.9, cex=1, sfrac=0, ylim=c(-0.005, 0.005), xlim=c(0.5,3.5), main="Number of Casualties")
abline(h=0, col="black")
axis(side=1, at=1:3, labels= labels.casualties, cex.axis=0.8, padj=1)
text(1.3, results_right_second_intif_casualties$estimate[1], labels=paste0(round(results_right_second_intif_casualties$pct_change[1], 2),"%"), col="red", cex=0.8)
text(2.3, results_right_second_intif_casualties$estimate[2], labels=paste0(round(results_right_second_intif_casualties$pct_change[2], 2),"%"), col="red", cex=0.8)
text(3.3, results_right_second_intif_casualties$estimate[3], labels=paste0(round(results_right_second_intif_casualties$pct_change[3], 2),"%"), col="red", cex=0.8)
dev.off()

# Figure 8: Placebo: Usage of Right-Wing Phrases in Different Specifications of "Peak" Years

file_names <- list.files("Data/Placebo_in_time/")
results_placebo <- list()
for (i in 1:length(file_names)){
	results_placebo[[i]] = read.dta(paste0("Data/Placebo_in_time/", file_names[i]))
}
names(results_placebo) <- file_names

results_treatment_right <- results_control_right <- results_treatment_center <- results_control_center <- results_treatment_left <- results_control_left <- as.data.frame(matrix(NA, nrow=(length(results_placebo)-8), ncol=ncol(results_placebo[[1]])))
colnames(results_treatment_right) <- colnames(results_control_right) <- colnames(results_treatment_center) <- colnames(results_control_center) <- colnames(results_treatment_left) <- colnames(results_control_left) <- colnames(results_placebo[[1]])


for (j in 1:(length(results_placebo)-8)){
	results_treatment_right[j,] <- results_placebo[[j]][results_placebo[[j]]$parm=="1.treatmentplacebo#1.right",]
	results_control_right[j,] <- results_placebo[[j]][results_placebo[[j]]$parm=="1.right",]
	
	results_treatment_center[j,] <- results_placebo[[j]][results_placebo[[j]]$parm=="1.treatmentplacebo#1.center",]
	results_control_center[j,] <- results_placebo[[j]][results_placebo[[j]]$parm=="1.center",]
	
	results_treatment_left[j,] <- results_placebo[[j]][results_placebo[[j]]$parm=="1.treatmentplacebo#1.left",]
	results_control_left[j,] <- results_placebo[[j]][results_placebo[[j]]$parm=="1.left",]

}
results_treatment_right$peak_year <- results_control_right$peak_year <- results_treatment_center$peak_year <- results_control_center$peak_year <- results_treatment_left$peak_year <- results_control_left$peak_year <- (1986:2006)

israelis.killed <- read.csv("Data/Israeli_casualties.csv", as.is=T)

empty_row <- c(1986, rep(NA, ncol(israelis.killed)-1))

israelis.killed <- rbind(empty_row, israelis.killed)
israelis.killed <- israelis.killed[-c(22,23),]

pdf("Figures/placebo_6years_right.pdf", width=9, height=5)
par(mar=c(5,4,3,4))
plot(israelis.killed$Year, israelis.killed$Total.killed, type="l", main="", xlab="", ylab="", xaxt="n", yaxt="n", col=NA)
rect(xleft=1987, ybottom=0, xright=1991, ytop=max(israelis.killed$Total.killed, na.rm=T)+10, col="lightgray", border=NA)
rect(xleft=2000, ybottom=0, xright=2005, ytop=max(israelis.killed$Total.killed, na.rm=T)+10, col="lightgray", border=NA)
lines(israelis.killed$Year, israelis.killed$Total.killed)
text(1989, 50, labels="First Intifada")
text(2002.5, 30, labels="Second Intifada")
axis(side=4, at=seq(0,450,100))
mtext("Number of casualties", side=4, line=2.5)
par(new=T)
plotCI(x= results_treatment_right$estimate, ui= results_treatment_right$max95, li= results_treatment_right$min95, type="p", pt.bg= "red", col="red", pch=21, gap=0, xaxt="n", ylab="Coefficient: Change in frequency of right-wing phrases", xlab="Placebo Peak Year", cex.main=1, cex.axis=0.9, cex=1, sfrac=0, ylim=c(-.0000002,.0000002))
abline(h=0, col="black")
axis(side=1, at=1:21, labels=as.character(results_treatment_right$peak_year), cex.axis=0.6)
dev.off()

# Figure A4: Residuals over time

resids = read.csv("Data/residuals_66model.csv", as.is=T)
resids2 = resids
resids2$bigram_year = paste(resids2$bigram, resids2$year, sep="_")
resids2 = resids2[duplicated(resids2$bigram_year)==FALSE,]
bigram_lengths = aggregate(resids2, FUN=length, by=list(bigram=resids2$bigram))
bigram_lengths = bigram_lengths[bigram_lengths[,2] == max(bigram_lengths[,2]),]
resids3 = resids2[resids2$bigram %in% bigram_lengths$bigram,]


pdf(file="Figures/scatter_resids_66_overtime.pdf", width=6, height=5)
plot(resids3$year[resids3$year>1995], resids3$resids_66[resids3$year>1995], pch=16, col=rgb(red=0.2, green=0.2, blue=1, alpha=0.2), ylab="Residuals", xlab="Year", cex.lab=0.8, cex.axis=0.8)
abline(h=0)
dev.off()

# Figure A5:  Year-by-year correlation of the residuals

pdf(file="Figures/scatter_resids_66_by_years.pdf", width=6, height=7)
par(mfrow=c(4,3))
par(mar=c(4,4,2,2))
plot(resids3$resids_66[resids3$year==1997], resids3$resids_66[resids3$year==1996], pch=16, col=rgb(red=0.2, green=0.2, blue=1, alpha=0.2), ylab="Residuals (1997)", xlab="Residuals (1996)", cex.lab=0.7, cex.axis=0.7, xlim=c(-400,400), ylim=c(-400,400))
par(mar=c(4,4,2,2))
plot(resids3$resids_66[resids3$year==1998], resids3$resids_66[resids3$year==1997], pch=16, col=rgb(red=0.2, green=0.2, blue=1, alpha=0.2), ylab="Residuals (1998)", xlab="Residuals (1997)", cex.lab=0.7, cex.axis=0.7, xlim=c(-400,400), ylim=c(-400,400))
par(mar=c(4,4,2,2))
plot(resids3$resids_66[resids3$year==1999], resids3$resids_66[resids3$year==1998], pch=16, col=rgb(red=0.2, green=0.2, blue=1, alpha=0.2), ylab="Residuals (1999)", xlab="Residuals (1998)", cex.lab=0.7, cex.axis=0.7, xlim=c(-400,400), ylim=c(-400,400))
par(mar=c(4,4,2,2))
plot(resids3$resids_66[resids3$year==2000], resids3$resids_66[resids3$year==1999], pch=16, col=rgb(red=0.2, green=0.2, blue=1, alpha=0.2), ylab="Residuals (2000)", xlab="Residuals (1999)", cex.lab=0.7, cex.axis=0.7, xlim=c(-400,400), ylim=c(-400,400))
par(mar=c(4,4,2,2))
plot(resids3$resids_66[resids3$year==2001], resids3$resids_66[resids3$year==2000], pch=16, col=rgb(red=0.2, green=0.2, blue=1, alpha=0.2), ylab="Residuals (2001)", xlab="Residuals (2000)", cex.lab=0.7, cex.axis=0.7, xlim=c(-400,400), ylim=c(-400,400))
par(mar=c(4,4,2,2))
plot(resids3$resids_66[resids3$year==2004], resids3$resids_66[resids3$year==2003], pch=16, col=rgb(red=0.2, green=0.2, blue=1, alpha=0.2), ylab="Residuals (2004)", xlab="Residuals (2003)", cex.lab=0.7, cex.axis=0.7, xlim=c(-400,400), ylim=c(-400,400))
par(mar=c(4,4,2,2))
plot(resids3$resids_66[resids3$year==2005], resids3$resids_66[resids3$year==2004], pch=16, col=rgb(red=0.2, green=0.2, blue=1, alpha=0.2), ylab="Residuals (2005)", xlab="Residuals (2004)", cex.lab=0.7, cex.axis=0.7, xlim=c(-400,400), ylim=c(-400,400))
par(mar=c(4,4,2,2))
plot(resids3$resids_66[resids3$year==2006], resids3$resids_66[resids3$year==2005], pch=16, col=rgb(red=0.2, green=0.2, blue=1, alpha=0.2), ylab="Residuals (2006)", xlab="Residuals (2005)", cex.lab=0.7, cex.axis=0.7, xlim=c(-400,400), ylim=c(-400,400))
par(mar=c(4,4,2,2))
plot(resids3$resids_66[resids3$year==2007], resids3$resids_66[resids3$year==2006], pch=16, col=rgb(red=0.2, green=0.2, blue=1, alpha=0.2), ylab="Residuals (2007)", xlab="Residuals (2006)", cex.lab=0.7, cex.axis=0.7, xlim=c(-400,400), ylim=c(-400,400))
par(mar=c(4,4,2,2))
plot(resids3$resids_66[resids3$year==2008], resids3$resids_66[resids3$year==2007], pch=16, col=rgb(red=0.2, green=0.2, blue=1, alpha=0.2), ylab="Residuals (2008)", xlab="Residuals (2007)", cex.lab=0.7, cex.axis=0.7, xlim=c(-400,400), ylim=c(-400,400))
dev.off()

# Figure A6: Placebo: Terrorism casualties from other countries

coeffs.results_gtd <- read.csv("Data/coeffs-for-graph-gtd.csv")

pdf(file="Figures/placebo_results_gtd.pdf", width=5, height=5)
plot(seq(1,6,1), coeffs.results_gtd$estimate, type = "p", axes = F, xlab = "", ylab = "", pch = 19, cex = 0.8, xaxs = "r", col="darkred", ylim=c(-0.0000000008, 0.0000000008))        
segments(seq(1,6,1), coeffs.results_gtd$X95min,seq(1,6,1), coeffs.results_gtd$X95max, lwd =  1, col="darkred")
axis(2, at = seq(-0.0000000008,0.0000000008,by= 0.0000000004), labels =seq(-0.0000000008,0.0000000008,by= 0.0000000004), tick = T,  cex.axis = 0.8, mgp = c(2,.7,0))  
axis(1, at = seq(1,6,1), labels =NA, tick = T,  cex.axis = 0.8, mgp = c(2,.7,0))
mtext(text= coeffs.results_gtd$X, las=1, side=1, at= seq(1,6,1), line=0.5, cex=0.7, padj=1)
abline(h=0)
title(main = "Placebo: Terrorism casualties from other countries", cex.main=1)
mtext("Change in phrase frequency", side=2, line=2.5, cex=0.8)
mtext("5-year moving average of casualties", side=1, line=2.5, cex=0.8)
dev.off()
